Scaling of Circulation in Buoyancy Generated Vortices 
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The temporal evolution of the fluid circulation generated by a buoyancy force when two- 
dimensional (2D) arrays of 2D thermals are released into a quiescent incompressible fluid is studied 
through the results of numerous lattice Boltzmann simulations. It is observed that the circulation 
magnitude grows to a maximum value in a flnite time. When both the maximum circulation and the 
time at which it occurs are non-dimensionalised by appropriately deflned characteristic scales, it is 
shown that two simple Prandtl number (Pr) dependent scaling relations can be devised that flt these 
data very well over nine decades of Pr spanning the viscous and diffusive regimes and six decades 
of Rayleigh number (Ra) in the low Ra regime. Also, obtained analytically is the exact result that 
circulation magnitude continues to grow in time for a two-dimensional laminar or turbulent single 
buoyant (3D) vortex ring in an inflnite unbounded fluid. 



Buoyancy generated vorticity is an attractive area of 
study for its influential role in various fields of science 
and engineering, its relevance to mixing and, undoubt- 
edly, for the aesthetically pleasingnature of its visualised 
flow structure (see, for example, [i|, [1,0, II, S 01 and the 
references therein) . Buoyant vortices are generated when 
plumes and thermals have different temperatures to an 
ambient fluid. Following some classical work on buoyant 
vortex rings and starting plumes 0,131, detailed studies 
on a single starting plume have generated some under- 
standing, though not yet complete, of the interesting phe- 
nomena of mushroom-type vortex head generation and 
its pinch-off (see, for example, Q). A number of studies 
on thermals have focused mainly on the time evolution of 
the linear dimension of the thermal and its penetration in 
the streamwise direction [l^, llll Il2| • Lundgren et al. [lfl| 
has also provided some information on the time evolution 
of the circulation. Further, effects of the initial geometry 
of thermals on their evolution were investigated recently 
[ij, |lj|. Despite this considerable interest, some basic 
aspects of buoyancy generated vorticity remain to be elu- 
cidated, including the presence or otherwise of quantita- 
tive scaling laws in different flow regimes. Such scaling 
laws may be used to predict aspects of the behaviour of a 
system without performing full solutions of the system's 
governing equations. 

In this letter, the results of numerous computer sim- 
ulations using the lattice Boltzmann method (LBM) are 
used to investigate the universal scaling behaviour asso- 
ciated with the circulation generated by buoyant forces 
when 2D thermals are released into a quiescent incom- 
pressible fluid. Also presented is an analytical derivation 
showing that for a 3D buoyant vortex ring formed by 
releasing a thermal in an inflnite domain of a quiescent 
fluid the magnitude of the circulation grows continuously 
in time for both laminar and turbulent cases. 

In terms of the Cartesian coordinate system {x, y, z) 
each of the simulated 2D systems comprised a 2Lx x 
2Ly X 1 sized domain of incompressible fluid which was 



initially quiescent and of uniform density, p, and tem- 
perature. To- The centre of the domain coincided with 
the origin of a Cartesian coordinate system. At time 
t — Q & circular (to within the lattice resolution) thermal 
of initial radius i?o < L^, Ly and temperature Ti < Tq 
was introduced into the centre of the domain. CycHc 
boundary conditions were applied at all lattice bound- 
aries making the simulation equivalent to that of an infl- 
nite system initialised with an inflnite number of circular 
thermals positioned at the nodes of a rectangular array 
such that their centres were separated by 2Lx in the x- 
direction and 2Ly in the y-direction. The temperature 
difference causes the thermals to move in the negative 
(downward) y-direction which coincides with the direc- 
tion of acceleration due to gravity —gey, where ey is the 
unit vector in the y-direction. This motion is caused by 
the buoyancy force acting on the thermal due to its den- 
sity being different from that of the surrounding fluid. 
For i > the temperature of the thermal diffuses and 
convects as it descends and a vorticity fleld with non- 
zero component uJz{x, y, t) = uj{x, y,t) — V A u{x, y, t) = 
dxUy{x,y,t) — dyUx(x,y,t) is generated, where Ux and 
Uy are the fluid velocity components in the x and y- 
directions, respectively. This system is governed by the 
Navier-Stokes equations with the Boussinesq approxima- 
tion, written as 



diUi = 0, 
DtUi = -p^^dip - SiyQ + ag{T - Tq) S., 

and the equation for temperature fleld T 

DtT = nd^T, 



ly 



vd'^Ui 



where summation notation applies to the indices i and j 
which take the values x and y, Dt = dt + Uidi, Sij is the 
Kronecker delta function, p is the pressure, v is the kine- 
matic viscosity and k is the thermal diffusivity. The solu- 
tion of these equations with cyclic boundary conditions 
at a; = Lx,~Lx and y — Ly,—Ly and initial tempera- 
tures Ti and Tq describes the flow fleld generated by the 



infinite rectangular array of thermals. It should be noted 
that Lx ^ oo and Ly ^ oo represents the situation of a 
single isolated thermal in an infinite fiuid. 

The LBM used to solve the governing equations was 
a multi-relaxation-time algorithm which sets all non- 
hydrodynamic modes to zero at each time step. It is 
based on the LBM for the Boussinesq equations de- 
scribed in |13]. At regular intervals during the simula- 
tions measurements were made of the circulation, T{t), 
calculated over the a; > half-domain and defined by 
r(i) = J_2 /g "^ ijj{x,y,t)dxdy. It was observed in all 
simulations that the circulation magnitude, r(i), grew 
to a maximum value, denoted Fniax, in a time denoted 
by ^max, and then decayed. The parameters infiuencing 
the phenomena are gaAT, v^ k, Rq, L^ and Ly (here, 
AT = \To — Ti\). In general 



i max 
imax 



^maxigaAT,u,K,Ro,Lx,Ly), (1) 

tmax{gaAT,iy,K,Ro,Lx,Ly). (2) 



To reveal the scaling relations for Fniax and imax a 
total of 211 individual simulations were performed us- 
ing various combinations of the parameters gaAT, v, k, 
Rq, L^ and Ly. Sets of simulations were performed in 
which all variables except one were held constant and 
the values Fniax and tmax were measured and recorded. 
The simulations were run for sufficient times, between 
3, 000 and 160, 000 time steps depending on the pa- 
rameter set, to enable the measurement of the initial 
growth of circulation and its subsequent decay. The pa- 
rameter ranges, in lattice Boltzmann units, were as fol- 
lows: 10-5 < agAT < 5 x 10-^ IQ-" < j/ < 19/6, 
10-"* < K < 1, 2/V7F < Ro < 4^5/^, 30 < L^ < 500 
and 30 < Ly < 500 ^^. These dimensional parameters 
where combined to yield dimensionless parameters cov- 
ering a wide range of values. Specifically: the Prandtl 
number, Pr = i^/k, ranged over nine decades (10^* < 
Pr < 10"*); the Rayleigh number, Ra = gaATR^/vn, 
over six decades (1.4 x 10"^ < Ra < 38.5); and the as- 
pect ratios A^^ = L^/Rq and Aj, = Ly/Ro, over nearly 
two decades (7.9 < Ax.y < 443). The range of Pr crosses 
from the low to high (diffusive to viscous) Pr regimes; 
the range of Ra, however, remains in the low Ra regime. 
Firstly, the dependency of the parameters gaAT, Rq, 
Lx or Ly was investigated by varying just one of these pa- 
rameters and examining plots of Fniax and imax versus 
the varied parameter. This showed the following scaling 
relations to hold: Fniax ^ gaAT ^ R^ ^ Lx ^ Ly and 

tmax ~ {gaATf r^ R^ ^ Ll r^ L". The dependency 
on i^ or K was more complicated, as varying just one 
of these and examining plots of Fniax and imax versus 
the varied parameter showed behaviour that depended 
on Pr. These relations suggest appropriate characteristic 
scales for circulation and time are Fo = LxRq gaAT/v 
and to = Lx/i^, where the viscosity is used in the denom- 
inator to ensure the correct dimensionality ( note that 
one could have equally well used k instead of i^). These 
characteristic scales, Fq and to, thus contain the correct 



dependency of Fniax and tmax on all the parameters ex- 
cept I' and K. Along with Q and (|2l, this suggests that 
the maximum circulation and the time of its occurrence 
may be written in the following dimensionless form which 
depends only on Pr 



Fmax/Fo 
tmax /to 



/i (Pr) 
/2 (Pr) 



(3) 
(4) 



The functions /i (Pr) and /2 (Pr) are as yet unknown, 
but, to see the scaling relations, Fniax/ro and tniax/to 
are plotted against Pr on log-log scales in Figure Q] The 
abscissa covers nine decades of Pr and the ordinate, 
about four decades of Fmax/Fo and tmax/to- The fig- 
ures exhibit certain power laws which may be written as 
Fmax/Fo c>c Pr" and tmax/to oc Pr"*, where the values 
of n and m are observed to be different in different Pr 
regimes. The plots in Figure Q] suggest that the values 
of n and m tend to constant values in the limits of high 
and low Pr. Assuming this to be the case, the following 
scaling relations, chosen for their simple form and correct 
asymptotic behaviour, were fitted to the data 



'max/Fo - 


= a(Pr-"'+Pr-"") \ 


(5) 


tmax/to = 


= 6(Pr-™'+Pr-"")"\ 


(6) 



where a, m, Uh, b, mi and mh are constant fitting param- 
eters. A nonhnear least-squares Marquardt-Levenberg 
algorithm was used to fit the data. The fitted scaling 
relations are shown as black lines on the plots in Fig- 
ure m and it can be seen that they do indeed give ex- 
tremely good fits to the data across all the Pr regimes. 
The fitted values of the fitting parameters are, with stan- 
dard errors: a = 0.4936 ± 0.0046, m = 0.9458 ± 0.0015, 
Uh = 0.0539 ± 0.0014, b = 0.20912 ± 0.0015, mi = 
0.8770 ± 0.0013 and mh = 0.1275 ± 0.0012. To within 
three standard errors all of these may be given as the fol- 
lowing simple fractions: a — 1/2, n; — 19/20, n^ = 1/20, 
b — 1 ji^, mi = 7/8 and mt — 1/8. It is an in- 
teresting observation, for which we have no explana- 
tion, that the fitted values suggest that n; = I — n^ 
and 771/ = 1 — m,ii. This suggests single exponent scal- 



Pr 



2n,, 



and 



ings of the form Fmax/Fo = a Pr'"- (1 

tmax/to = 6 Pr"''(l + Pr2'"''-i)"\ 

The dependencies on the other dimensionless parame- 
ters, Ra and the aspect ratio A^; = Lx/Ro, can be high- 
lighted by writing the scalings relations IJSI and 1^ as 
Fmax/K = A^Ra/i(Pr) and i^tmax/Ro = A2/2(Pr). 
The accuracy of scalings obtained for Fmax and tmax 
are further confirmed by the plots in Figures El and El 
of the alternatively non-dimensionalised Fmax and tmax 
plotted against the two dimensionless parameters Ra and 
Ax- In both figures the fitted scaling relations are shown 
as black lines, and again, the fits are seen to be extremely 
good. 

The scaling relations given in ij^l and (JSJ suggest that 
in an infinite domain, that is to say, Lx —>■ oo, Fmax 
and tmax tend to infinity. A theoretical justification for 
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Figure 1: (Colour online) Dimensionless Fmax/ro and 
imax/io plotted over six decades of Prandtl number. The 
black lines are the best fits to the functions using Q and ©• 



this observation is now provided. The equation for the 
non-zero component of vorticity, lu, can be written as 

dtUJ + dx {uxUj) + dy {uyUj) = v [d^ + 9y) u + agdxT 

and so the equation governing the circulation in an infi- 
nite domain can be written as 

dtV = ag / (To - T\x=o) dy - v I {dxu:\x=Q) dy. 



(7) 
The first term on the right-hand side (rhs) of 10 is the 
buoyancy term, which is always positive as Tq > Ti forces 
(To — T\x=o) > 0. The second term on the rhs of 10 is 
negative due to the fact that dxOj\x=o > because of the 
change in the sign of uj near a: = 0. Also, the magni- 
tude of this second term is zero at t = and starts to 
increase as the vorticity field is generated in the domain. 
In the limit of the viscous force being small compared to 
the buoyancy force, dtT remains positive and so F will 
continue to increase without limit. Although the above 
argument for a continuous increase of F requires v to be 
small, this is not a restriction in the case of an axisym- 
metric buoyant vortex ring in an infinite domain as is 
now discussed. 
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Figure 2: (Colour online) Particular non-dimensionalised 
forms of Fmax and tmax plotted over six decades of Rayleigh 
number. The black lines are Fmax/«Aa;/i (Pr) = Ra and 
i^imax/A^i?o.f2 (Pr) = 1 respectively. 



Consider a domain of quiescent fiuid of uniform density 
p and temperature Tq in a cylindrical coordinate system 
(r, 6, z) having < r < cx), — oo < 2; < oo and < 6* < 27r 
and acceleration due to gravity, 5, acting in the negative 
z-direction. At time i = a spherical volume of fiuid of 
radius Rq and temperature Ti < Tq is introduced into the 
domain with its centre coinciding with the origin of the 
coordinate system. For time i > 0, the thermal starts 
descending along the z-axis due to the buoyancy force 
and a vortex ring of buoyant thermal fiuid is generated 
(see, for example, |3| for a discussion of an ascending 
vortex ring when Ti > To). Due to the angular sym- 
metry in the 6'-direction, this system can be analysed in 
2D (r, z) coordinates. By employing the Boussinesq ap- 
proximation for the buoyancy force in the Navier-Stokes 
equations and using uje — dzUr — drUz to denote the non- 
zero vorticity component, the governing equation for the 
circulation T{t) — J^ J_^ LUg dz dr of the 2D vortex ring 
can be written as 



dtT ^ ag {T{r = 0, z, t) - Tq) dz, 



(8) 



in which the viscous term does not arise due to the sym- 
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Figure 3: (Colour online) Particular non-dimensionalised 
forms of Fmax and tmax plotted over two decades of the 
aspect ratio A^,. The black lines are Fmax/KRa/i (Pr) — A^ 
and f tmax/i?o h (Pr) = ^l respectively. 



dt{r} = agJ'^^{{T{r^O,z,t)}-To)dz, where (F) = 

la '^^ /_oo '^'^ ^^^) ^^ ^^^ ensemble average of the instan- 
taneous circulation F = (F) + F' and F' represents the 
turbulent fluctuations in F over the (F). Thus, because 
T{r = 0, z, t)-To < 0, or {T{r = 0, z, t)) - Tq < for a 
turbulent system, (jHJ suggests that in the case of a de- 
scending vortex ring in an infinite domain both F and 
(F) continue to decrease (or |r| and |(r)| continue to in- 
crease) from their initial values of zero. 

In this letter it was shown that, for an infinite sys- 
tem comprising an infinite number of 2D thermals ini- 
tially arranged in a rectangular array, the magnitude of 
the buoyancy generated circulation, F(i), reaches a max- 
imum value, Fmax, at a finite time, imax- Accurate scal- 
ing relations for Fmax and imax, covering nine decades 
of Pr and six decades of Ra, were inferred from LB sim- 
ulations. Theoretical justification has been provided to 
support the observation, based on the scaling relations, 
that Fmax increases in proportion to the size of the do- 
main. Furthermore, exact analytical results were derived 
for a single buoyancy generated vortex ring in both the 
laminar and turbulent cases. These exact results suggest 
that for a buoyant vortex ring in an infinite unbounded 
domain the magnitude of F, or (F) if the vortex ring is 
turbulent, will continue to grow indefinitely. The impli- 
cation of this in predicting the growth of size R of the 
vortex ring can be exhibited by the second term involv- 
ing dT/dt in npT dR^/dt + npR^dT/dt = Fb, where Fb is 
constant buoyancy force acting on the buoyant fiuid. In 
fact. Turner Q assumed F was constant and neglected 
this second term when predicting the growth of a vortex 
We also note that it is assumed in the model used in 
that the circulation remains constant after an initial 



metry of the problem. Indeed, starting from the en- 
semble averaged (denoted by (•)) Navier-Stokes equa- 
tions and applying various symmetry conditions to the 
ensemble averaged fiow properties, one can show that 



It is hoped that these results will inspire future work, 
for example: finding scaling relations and exponents in 
3D cases with different initial configurations of thermals, 
addressing why the exponents found here add up to unity, 
and extending the study to examine high Ra regimes. 
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